#ifndef _TESTMULTIGRID_H_
#define _TESTMULTIGRID_H_

#include "Multigrid/MGFactory.h"
#include <ctime>

template <int Dim, class RST, class ITPL> class TestMultigrid;

template <int Dim, class RST, class ITPL>
class TestMultigrid{
public:

  using iVec = Vec<int,2>;

  using rVec = Vec<Real,2>;
  
  template <class T>
  using Vector = std::vector<T>;

  TestMultigrid(const MGFactory<Dim,RST,ITPL>& aMGF);

  template <class TFunc>
  Vector<Real> encrypTestVCycle(const TFunc& f, const int numEncryp, const
                          int p = 0);

  template<class TFunc>
  void TestVCycle(const TFunc &f, const int numEncryp);

  template<class TFunc>
  void TestVCycleCPUTimes(const TFunc &f, const int numEncryp);

  template<class TFunc>
  void TestFMVCycle(const TFunc &f, const int numEncryp);
  
protected:
  MGFactory<Dim,RST,ITPL> MGF;
};

template <int Dim, class RST, class ITPL>
TestMultigrid<Dim,RST,ITPL>::TestMultigrid(const MGFactory<Dim,RST,ITPL>& aMGF):MGF(aMGF){}

template <int Dim, class RST, class ITPL>
template<class TFunc>
std::vector<Real> TestMultigrid<Dim,RST,ITPL>::encrypTestVCycle(const
                                                                TFunc &f, const int numEncryp, const int p){
  std::vector<Real> res(2*numEncryp+1);
  const char* BCTypes = MGF.getBCTypes();
  const RectDomain<Dim>& Domain = MGF.getDomain();
  int n = Domain.size()[0]-1;
  const Real weight = MGF.getWeight();
  const Vec<int,3> numIter = MGF.getNumIter();
  const Real reltol = MGF.getReltol();
  const int maxIter = MGF.getMaxIter();
  switch(Dim){
  case 1:
    std::cout << "Grid " << n << ":" << std::endl;
    break;
  case 2:
    std::cout << "Grid " << n << "X" << n << ":" << std::endl;
    break;
  }
  res[0] = MGF.computeError(MGF.SolveVCycle(f),f,p);
  std::vector<RectDomain<Dim> > vDomain{Domain};
  for(int i = 1 ; i <= numEncryp ; i++){
    RectDomain<Dim> nDomain = vDomain.back().refine();
    vDomain.push_back(nDomain);
    MGFactory<Dim,RST,ITPL> nMGF(nDomain,BCTypes);
    nMGF.setParam(weight,numIter[0],numIter[1],numIter[2],reltol,maxIter);
    n*=2;
    switch(Dim){
    case 1:
      std::cout << "Grid " << n << ":" << std::endl;
      break;
    case 2:
      std::cout << "Grid " << n << "X" << n << ":" << std::endl;
      break;
    }
    res[2*i] = nMGF.computeError(nMGF.SolveVCycle(f),f,p);
    res[2*i-1] = log(res[2*i-2]/res[2*i])/log(2);
  }
  return res;
}

template <int Dim, class RST, class ITPL>
template<class TFunc>
void TestMultigrid<Dim,RST,ITPL>::TestVCycle(const TFunc &f, const int numEncryp){
  std::vector<Real> res1(2*numEncryp+1);
  std::vector<Real> res2(2*numEncryp+1);
  std::vector<Real> res3(2*numEncryp+1);
  const char* BCTypes = MGF.getBCTypes();
  const RectDomain<Dim>& Domain = MGF.getDomain();
  int n = Domain.size()[0]-1;
  const Real weight = MGF.getWeight();
  const Vec<int,3> numIter = MGF.getNumIter();
  const Real reltol = MGF.getReltol();
  const int maxIter = MGF.getMaxIter();
  switch(Dim){
    case 1:
      std::cout << "Grid " << n << ":" << std::endl;
      break;
    case 2:
      std::cout << "Grid " << n << "X" << n << ":" << std::endl;
      break;
  }
  Tensor<Real,Dim> phi = MGF.SolveVCycle(f);
  res1[0] = MGF.computeError(phi,f,0);
  res2[0] = MGF.computeError(phi,f,1);
  res3[0] = MGF.computeError(phi,f,2);
  std::vector<RectDomain<Dim> > vDomain{Domain};
  for(int i = 1 ; i <= numEncryp ; i++){
    RectDomain<Dim> nDomain = vDomain.back().refine();
    vDomain.push_back(nDomain);
    MGFactory<Dim,RST,ITPL> nMGF(nDomain,BCTypes);
    nMGF.setParam(weight,numIter[0],numIter[1],numIter[2],reltol,maxIter);
    n*=2;
    switch(Dim){
    case 1:
      std::cout << "Grid " << n << ":" << std::endl;
      break;
    case 2:
      std::cout << "Grid " << n << "X" << n << ":" << std::endl;
      break;
    }
    Tensor<Real,Dim> cphi = nMGF.SolveVCycle(f);
    res1[2*i] = nMGF.computeError(cphi,f,0);
    res2[2*i] = nMGF.computeError(cphi,f,1);
    res3[2*i] = nMGF.computeError(cphi,f,2);
    res1[2*i-1] = log(res1[2*i-2]/res1[2*i])/log(2);
    res2[2*i-1] = log(res2[2*i-2]/res2[2*i])/log(2);
    res3[2*i-1] = log(res3[2*i-2]/res3[2*i])/log(2);
  }
  std::cout << "----------------------------------------" << std::endl;
  const Vec<Real,Dim>& dx = Domain.spacing();
  const Vec<int,Dim>& lo = Domain.lo();
  const Vec<int,Dim>& hi = Domain.hi();
  n = hi[0] - lo[0];
  switch(Dim){
  case 1:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")$" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] <<
    std::endl;
    break;
  case 2:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")\\times(" << lo[1] << "," << (hi[1]-lo[1])*dx[1] <<  \
      ")" << "$\\\\" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] << " , " \
              << BCTypes[2] << " , " << BCTypes[3] << std::endl;
    break;
  }
  std::cout << "\\begin{table}[htbp]" << std::endl;
  std::cout << "\\centering\\begin{tabular}{c|";
  for (int i = 1 ; i <= 2*numEncryp+1 ; i++)
    std::cout << "c";
  std::cout << "}" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$n$&" << n;
  for (int i = 1 ; i <= numEncryp ; i++){
    n*=2;
    std::cout << "&ratio&" << n;
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_1$&" << res2[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res2[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_2$&" << res3[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res3[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_{\\infty}$&" << res1[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res1[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "\\end{tabular}" << std::endl;
  std::cout << "\\caption{}" << std::endl;
  std::cout << "\\end{table}" << std::endl;

}

template <int Dim, class RST, class ITPL>
template<class TFunc>
void TestMultigrid<Dim,RST,ITPL>::TestFMVCycle(const TFunc &f, const int numEncryp){
  std::vector<Real> res1(2*numEncryp+1);
  std::vector<Real> res2(2*numEncryp+1);
  std::vector<Real> res3(2*numEncryp+1);
  const char* BCTypes = MGF.getBCTypes();
  const RectDomain<Dim>& Domain = MGF.getDomain();
  int n = Domain.size()[0]-1;
  const Real weight = MGF.getWeight();
  const Vec<int,3> numIter = MGF.getNumIter();
  const Real reltol = MGF.getReltol();
  const int maxIter = MGF.getMaxIter();
  // switch(Dim){
  //   case 1:
  //     std::cout << "Grid " << n << ":" << std::endl;
  //     break;
  //   case 2:
  //     std::cout << "Grid " << n << "X" << n << ":" << std::endl;
  //     break;
  // }
  Tensor<Real,Dim> phi = MGF.SolveFMVCycle(f);
  res1[0] = MGF.computeError(phi,f,0);
  res2[0] = MGF.computeError(phi,f,1);
  res3[0] = MGF.computeError(phi,f,2);
  std::vector<RectDomain<Dim> > vDomain{Domain};
  for(int i = 1 ; i <= numEncryp ; i++){
    RectDomain<Dim> nDomain = vDomain.back().refine();
    vDomain.push_back(nDomain);
    MGFactory<Dim,RST,ITPL> nMGF(nDomain,BCTypes);
    nMGF.setParam(weight,numIter[0],numIter[1],numIter[2],reltol,maxIter);
    n*=2;
    // switch(Dim){
    // case 1:
    //   std::cout << "Grid " << n << ":" << std::endl;
    //   break;
    // case 2:
    //   std::cout << "Grid " << n << "X" << n << ":" << std::endl;
    //   break;
    // }
    Tensor<Real,Dim> cphi = nMGF.SolveFMVCycle(f);
    res1[2*i] = nMGF.computeError(cphi,f,0);
    res2[2*i] = nMGF.computeError(cphi,f,1);
    res3[2*i] = nMGF.computeError(cphi,f,2);
    res1[2*i-1] = log(res1[2*i-2]/res1[2*i])/log(2);
    res2[2*i-1] = log(res2[2*i-2]/res2[2*i])/log(2);
    res3[2*i-1] = log(res3[2*i-2]/res3[2*i])/log(2);
  }
  std::cout << "----------------------------------------" << std::endl;
  const Vec<Real,Dim>& dx = Domain.spacing();
  const Vec<int,Dim>& lo = Domain.lo();
  const Vec<int,Dim>& hi = Domain.hi();
  n = hi[0] - lo[0];
  switch(Dim){
  case 1:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")$" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] <<
    std::endl;
    break;
  case 2:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")\\times(" << lo[1] << "," << (hi[1]-lo[1])*dx[1] <<  \
      ")" << "$\\\\" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] << " , " \
              << BCTypes[2] << " , " << BCTypes[3] << std::endl;
    break;
  }
  std::cout << "\\begin{table}[htbp]" << std::endl;
  std::cout << "\\centering\\begin{tabular}{c|";
  for (int i = 1 ; i <= 2*numEncryp+1 ; i++)
    std::cout << "c";
  std::cout << "}" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$n$&" << n;
  for (int i = 1 ; i <= numEncryp ; i++){
    n*=2;
    std::cout << "&ratio&" << n;
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_1$&" << res2[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res2[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_2$&" << res3[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res3[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_{\\infty}$&" << res1[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res1[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "\\end{tabular}" << std::endl;
  std::cout << "\\caption{}" << std::endl;
  std::cout << "\\end{table}" << std::endl;

}

template <int Dim, class RST, class ITPL>
template<class TFunc>
void TestMultigrid<Dim,RST,ITPL>::TestVCycleCPUTimes(const TFunc &f, const int numEncryp){
  std::vector<Real> res1(2*numEncryp+1);
  std::vector<Real> res2(2*numEncryp+1);
  std::vector<Real> res3(2*numEncryp+1);
  std::vector<Real> time(2*numEncryp+1);
  clock_t begin, end;
  const char* BCTypes = MGF.getBCTypes();
  const RectDomain<Dim>& Domain = MGF.getDomain();
  int n = Domain.size()[0]-1;
  const Real weight = MGF.getWeight();
  const Vec<int,3> numIter = MGF.getNumIter();
  const Real reltol = MGF.getReltol();
  const int maxIter = MGF.getMaxIter();
  switch(Dim){
    case 1:
      std::cout << "Grid " << n << ":" << std::endl;
      break;
    case 2:
      std::cout << "Grid " << n << "X" << n << ":" << std::endl;
      break;
  }
  begin = clock();
  Tensor<Real,Dim> phi = MGF.SolveVCycle(f);
  end = clock();
  res1[0] = MGF.computeError(phi,f,0);
  res2[0] = MGF.computeError(phi,f,1);
  res3[0] = MGF.computeError(phi,f,2);
  time[0] = (double)(end - begin) / CLOCKS_PER_SEC;
  std::vector<RectDomain<Dim> > vDomain{Domain};
  for(int i = 1 ; i <= numEncryp ; i++){
    RectDomain<Dim> nDomain = vDomain.back().refine();
    vDomain.push_back(nDomain);
    MGFactory<Dim,RST,ITPL> nMGF(nDomain,BCTypes);
    nMGF.setParam(weight,numIter[0],numIter[1],numIter[2],reltol,maxIter);
    n*=2;
    switch(Dim){
    case 1:
      std::cout << "Grid " << n << ":" << std::endl;
      break;
    case 2:
      std::cout << "Grid " << n << "X" << n << ":" << std::endl;
      break;
    }
    begin = clock();
    Tensor<Real,Dim> cphi = nMGF.SolveVCycle(f);
    end = clock();
    res1[2*i] = nMGF.computeError(cphi,f,0);
    res2[2*i] = nMGF.computeError(cphi,f,1);
    res3[2*i] = nMGF.computeError(cphi,f,2);
    time[2*i] = (double)(end - begin) / CLOCKS_PER_SEC;
    res1[2*i-1] = log(res1[2*i-2]/res1[2*i])/log(2);
    res2[2*i-1] = log(res2[2*i-2]/res2[2*i])/log(2);
    res3[2*i-1] = log(res3[2*i-2]/res3[2*i])/log(2);
    time[2*i-1] = -log(time[2*i-2]/time[2*i])/log(2);
  }
  std::cout << "----------------------------------------" << std::endl;
  const Vec<Real,Dim>& dx = Domain.spacing();
  const Vec<int,Dim>& lo = Domain.lo();
  const Vec<int,Dim>& hi = Domain.hi();
  n = hi[0] - lo[0];
  switch(Dim){
  case 1:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")$" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] <<
    std::endl;
    break;
  case 2:
    std::cout << "\\noindent Domain: " << "$(" << lo[0] << "," << (hi[0]-lo[0])*dx[0] \
              << ")\\times(" << lo[1] << "," << (hi[1]-lo[1])*dx[1] <<  \
      ")" << "$\\\\" << std::endl;
    std::cout << "BCtype : " << BCTypes[0] << " , "<< BCTypes[1] << " , " \
              << BCTypes[2] << " , " << BCTypes[3] << std::endl;
    break;
  }
  std::cout << "\\begin{table}[htbp]" << std::endl;
  std::cout << "\\centering\\begin{tabular}{c|";
  for (int i = 1 ; i <= 2*numEncryp+1 ; i++)
    std::cout << "c";
  std::cout << "}" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$n$&" << n;
  for (int i = 1 ; i <= numEncryp ; i++){
    n*=2;
    std::cout << "&ratio&" << n;
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_1$&" << res2[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res2[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_2$&" << res3[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res3[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "$\\|\\mathrm{E}\\|_{\\infty}$&" << res1[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << res1[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "CPU time(s)&" << time[0];
  for (int i = 1 ; i < 2*numEncryp+1 ; i++){
    std::cout << "&" << time[i];
  }
  std::cout << "\\\\" << std::endl;
  std::cout << "\\hline" << std::endl;
  std::cout << "\\end{tabular}" << std::endl;
  std::cout << "\\caption{}" << std::endl;
  std::cout << "\\end{table}" << std::endl;
}



#endif // _TESTMULTIGRID_H_
